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Abstract. Several modifications aimed at improving physical accuracy are proposed for solving axially symmetric 
problems building on the DSMC (2007) algorithms introduced by Bird. Originally developed to solve nonequilibrium, 
rarefied flows, the DSMC method is now regularly used to solve complex problems over a wide range of Knudsen 
numbers. These new algorithms include features such as nearest neighbor collisions excluding the previous collision 
partners, separate collision and sampling cells, automatically adaptive variable time steps, a modified no-time counter 
procedure for collisions, and discontinuous and event-driven physical processes. Axially symmetric solutions require 
radial weighting for the simulated molecules since the molecules near the axis represent fewer real molecules than those 
farther away from the axis due to the difference in volume of the cells. In the present methodology, these radial 
weighting factors are continuous, linear functions that vary with the radial position of each simulated molecule. It is 
shown that how one defines the number of tentative collisions greatly influences the mean collision time near the axis. 
The method by which the grid is treated for axially symmetric problems also plays an important role near the axis, 
especially for scalar pressure. A new method to treat how the molecules are traced through the grid is proposed to 
alleviate the decrease in scalar pressure at the axis near the surface. Also, a modification to the duplication buffer is 
proposed to vary the duplicated molecular velocities while retaining the molecular kinetic energy and axially symmetric 
nature of the problem. 


INTRODUCTION 

With the recent introduction of modified direct simulation Monte Carlo (DSMC) algorithms 1 , simulations can be 
made more quickly using fewer numbers of simulated particles for some cases. Also, since fewer numbers of 
particles are required for a given flow, the range of Knudsen numbers readily simulated can be extended further into 
the continuum regime. The new algorithms are another step towards improving DSMC algorithms focusing on 
physical accuracy while improving efficiency. Some improvements in the past have included additions such as 
conservative species weighting 2 and the use of virtual sub-cells 3 . However, physical improvements to the simulation 
of axially symmetric flows are not prevalent in the literature. One problem encountered in the calculation of the 
number of tentative collisions is how to define the center of the collision cell when performing collisions. If care is 
not taken, the number of collisions near the axis may be too high, resulting in a higher mean collision rate and lower 
mean free path near the axis. 

Another important issue when performing axially symmetric simulations is the inclusion of a duplication buffer. 
The duplication buffer is used to reduce the likelihood that when a molecule moves towards the axis and is 
duplicated, it does not collide with a clone of itself. The clone is therefore placed in the duplication buffer and a 
different molecule from the buffer is placed into the flow. A modification to the duplication buffer is proposed that 
uses a collision cell duplication buffer rather than a global duplication buffer. This is done to promote a time- 
accurate conservation of mass within the collision cells. Another modification to the cloned molecules is presented 
where the radial and axial velocities are unchanged while the out-of-plane velocity is reversed, thus preserving the 
kinetic energy of the molecule duplicated while changing nothing about the axially symmetric nature of the problem. 

It will also be shown that the method by which one treats the grid and molecular movements plays an important 
role in the calculation of pressure along the axis. Comparisons will be made of two traditional methods to move the 
molecules through their respective grid scheme, each having benefits and detractors. A new, hybrid grid /molecule 
translation scheme will be proposed combining the benefits of both traditional schemes. Results of the new scheme 
will be shown to have a much better distribution of pressure along the axis than one traditional method combined 



with a lower grid merit parameter (mean collision separation divided by mean free path ( mcs/mfp )) throughout the 
flow field compared to the other traditional method. 

The baseline computations are performed using the DS2 4 code and the DSMC Analysis Code (DAC) 5,6 . 
However, most of the computations presented herein are results from the author’s code. When a specific code is not 
mentioned in the text, the author is referring to his own code. 

TREATMENT OF AXIALLY SYMMETRIC PROBLEMS 


A complete discussion of the most common treatment of axially symmetric problems can be found in Ref. 7. 
Only the relevant points to this discussion will be given in overview. One problem associated with DSMC 
computations for axially symmetric flows is the very small number of molecules near the axis (small volume cells) 
versus the large number of molecules far away from the axis (large volume cells) if the weight of each molecule is 
held constant. This has led to the introduction of radial weighting factors such that a simulated molecule located far 
from the axis represents more real molecules than one near the axis. The most recent definition of the radial 
weighting factor, W F , is: 

W F = 1 + r ^WF /^max W 


where r is the radial position of the molecule, y max is the maximum radius of the flow-field, and R WF defines how the 
weighting factor varies throughout the flow field (recommended value of 10,000 4 ). This results in a more uniform 
number of molecules per cell (for equal area cells in the x-y plane), but there is now a probability that each molecule 
will have to be either duplicated or deleted after each movement. Because of the duplication of molecules through 
the use of radial weighting, there would be identical molecules present at the exact same location if a duplication 
buffer didn’t exist. Therefore, when a molecule moves towards the axis and is duplicated, a copy of the molecule is 
placed in the duplication buffer and another molecule is taken from the buffer and inserted into the flow. This has 
traditionally been done using a global duplication buffer for the entire flow field. Weighting factors may be based 
either on the radius of the cell in which the molecule resides or on the radius of the molecule itself. It has been 
recommended 7 that cell-based factors not be used. Using the traditional axially symmetric methodologies described 
in Ref. 7, the main issues remaining for the new DSMC (2007) algorithms are the mean collision time (met) and 
pressure (P) along the axis. 

Throughout this paper, comparisons between different algorithms will be made using a common test case. A 10- 
cm radius sphere is used in a flow of non-reacting, molecular Nitrogen where both rotational and vibrational modes 
have been ignored. The free stream properties are a number density of 1.2e21 m" 3 , a free stream velocity of 5100 m/s 
and a free stream temperature of 190 K resulting in a Knudsen number of 0.01. The surface temperature was 500 K 
with the total number of simulated molecules around 600,000. It should be noted that none of the solutions 
presented are grid converged within the shock layer in the sense that the maximum value of the merit parameter, in 
each collision cell is not on the order of 0.1 or less throughout the entire flow field. However, for gross flow field 
features such as local mean collision time and scalar pressure, grid 
convergence should not be required for reasonable results. Figure 1 
presents the surface pressure computed using DS2 showing a 
pronounced decrease in pressure at the stagnation point (0 = 0-deg). 

DS2 implements all of the features of the modified algorithms 
presented in Ref. 1. Figure 2 presents the flow field P, met and 
mcs/mfp computed using DS2. This figure exemplifies the issues 
discussed above. The decrease in scalar pressure along the axis can 
be seen in Fig. 2. a. The decrease in mean collision time can be seen 
in Fig. 2.b. Figure 2.c presents the mcs/mfp. The maximum value 
of met/mfp for this simulation is around 0.8 at the surface. The 
average value in the free-stream is on the order of 0.1 -0.2. It is the 
purpose of the remainder of this paper to propose and examine 
methods to remove these non-physical behaviors for axially 
symmetric flows produced by the DSMC method while maintaining 
a relatively low merit parameter. 
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FIGURE 1 . Surface pressure from DS2. 
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FIGURE 2. Solution computed by DS2. 


CALCULATION OF NUMBER OF TENTATIVE COLLISIONS 


The number of tentative collisions ( Nq ) between similar molecules is calculated as: 

N c =y 2 n(N-l)(og) max At (2) 

where n is the number density, N is the number of molecules in the collision cell, (crg) max is the maximum value of 
the product of collision cross-section (a) and relative velocity (g), and At is the local time step. The main issue is in 
the definition of the number density in the collision cell. The number density can be defined as: 

n = NW F F N V~ 1 (3) 

where F N is the global ratio of real to simulated molecules, V c is the volume of the collision cell, and W F is an 
average value of the radial weighting factor where all radial weights are assumed to be an average value for the 
collision cell. It is the definition of the center of the collision cell that can be rather ambiguous for axially 
symmetric simulations. There have been two traditional ways to compute the center of the collision cell: 
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where r max is the maximum radial coordinate of the collision cell and r min is the minimum radial coordinate of the 
collision cell. As it turns out, both of these definitions result in the mean collision time being too small along the 
axis, as shown in Fig. 3. A more accurate definition of the number density at each time step can be written: 
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FIGURE 3. Computed mean collision times using three methodologies of computing number density within a collision cell. 


where no assumption has been made as to the weights of the molecules. The met using eq. (6) is presented in Fig. 
3.c. However, there still appears to be an error along the axis where the met is still slightly too low. The mean 
collision time within a sample cell can be calculated as: 


met = | 
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where N mo i s is the total number of sampled molecules, AT samp is the total sample time, N co u is the total number of 
collisions in the sample and N samp is the total number of samples. In equation 6, no assumption was made as to the 
weights of the molecules. In equation 7, the number of molecules and number of collisions are assumed to have a 
constant weight within the sample cell. Therefore, it depends on how one counts the number of molecules as well as 
the number of collisions as to what the mean collision time is calculated to be. It is also believed that the 
assumption that the simulated molecules are of the same weight when collisions are performed contributes to the 
issues discussed herein. It is the focus of future research to include the variation of molecular weights in the 
collision process. For now, the mean collision time has been improved, thus the original problem has been 
alleviated. The result of the modified method of calculating the number of tentative collisions can be seen in Fig. 4. 
Although the met has been improved, the scalar pressure along the axis remains low. 




FIGURE 4. Solution using modified number of tentative collisions. 


GRID SYSTEM AND MOLECULE TRANSLATION 

There are two methods discussed presently to treat axially symmetric flows in DSMC. The first method, 
employed by DS2, is where the molecules are moved in a three-dimensional manner and then rotated back to the 
symmetry plane, accounting for the rotation in both the position and velocity 7 . The other method is to treat the grid 
as a three-dimensional wedge centered about the symmetry plane where the ±z-faces of the wedge are treated as 
specular reflection surfaces, thus enforcing the symmetry. However, the molecular position and velocity are never 
rotated back to the symmetry plane, but the y- and z-velocities are changed when a reflection occurs. DAC treats 
axially symmetric problems in this manner with a 5 -degree wedge about the symmetry plane and it is based on 
Bird’s 1994 algorithms 7 , is not time accurate for grid adapted problems, and the radial weights of molecules are cell- 
based, not radial position-based. For the modeled conditions, pressure is rather 
well behaved near the axis for a grid- converged solution as shown in Fig. 5. 

This grid system was has been implemented using Bird’s 2007 algorithms 1 , 
now in a time-accurate implementation where the radial weighting factors vary 
continuously as a function of molecular position. The pressure field results of 
this implementation can be seen in Fig. 6. a and is markedly better along the 
axis than that in Fig. 4. a. However, the quality of the solution decreases as the 
radial position increases, especially along the surface. This is because of 
decreased resolution caused by the three dimensional nature of the grid. As the 
radial position increases, the width of the cells also increases. The values of 
mes/mfp are presented for this case in Fig. 6.c. Clearly, this parameter is 
increasing radially, resulting in a less accurate solution. 

It may be that part of the improvement in pressure distribution is because 



FIGURE 5. DAC pressure distribution. 




FIGURE 6. Results of simulation using a 5 -degree wedge about the symmetry plane, 
the specular reflections of the molecules reverse the component of the velocity vector perpendicular to the ±2.5- 
degree planes. Near the axis, where there may be multiple duplicates of molecules because of the radial weighting 
factor, this plays an important role in diversifying the population of molecules. Therefore, two changes are 
proposed. First, a modification to the molecule duplication process is proposed. When a molecule is duplicated, the 
velocity in the plane perpendicular to the radial direction and x-axis is multiplied by -1.0” c , where nc is the cloned 
molecule number for multiply duplicated molecules. While this does nothing to the axially symmetric nature of the 
problem, it adds diversity in the population. The pressure is defined as the product of the number density, 
Boltzmann’s constant and the translational temperature. Translational temperature is defined as: 
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Reversing the direction of the out-of-plane velocity will tend to decrease the value of W , which will in turn 
increase the translational temperature and thus the pressure. 

Secondly, another treatment of axially symmetric grid system/molecule translation is proposed to take advantage 
of the improvements gained by the specularly reflected treatment of molecular movement, such as with the DAC 
implementation. In order to decrease and thus improve thejnerit parameter, the angle of the reflection planes is 
transformed as a function of collision cell radial location. The reflection plane angle is defined as: 
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where r max and r min are the maximum and minimum radial coordinates of the collision cell, such that the maximum 
width of the collision cell in the z-direction is equivalent to the collision cell height at the maximum angle. This 
gives the grid a more three dimensional appearance while limiting the width of the cells. At the end of each 
molecular move, since the position in the angular-direction about the axis does not affect the axially symmetric 
solution, the position and velocity of the molecules are rotated to a random angular location within the collision cell 
it resides in (the maximum and minimum defined as ±|3). This 
ensures that the molecule is within the collision cell width and there 
is diversity in the molecule population without affecting the 
physical nature of the solution. The surface pressure computed by 
this method is presented in Fig. 7 with a comparison to the DS2 
results presented earlier. Using the modifications presented herein, 
there is still a slight decrease in pressure at the stagnation point, but 
not as pronounced. As stated earlier, it is believed that if the 
variation in molecular weights was accounted for in the collision 
process, this decrease would disappear. The flow field results are 
presented in Fig. 8. The pressure distribution along the axis closely 
resembles that of Fig. 5, and the mcs/mfp , while being larger than 
that of Fig. 4 (maximum of 1.49 versus 0.77), is much less than that 

of Fig. 6 (maximum of 1.49 versus 4.88). Angle (rad) 

FIGURE 7. Surface pressure comparison. 
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FIGURE 8. Solution using modified number of tentative collisions and reversal of out-of-plane velocity in duplicated 
molecules with specular reflection off of ±|3 planes off symmetry plane. 


CONCLUSION 

With the recent introduction of modified direct simulation Monte Carlo (DSMC) algorithms 1 , simulations can be 
made more quickly using fewer numbers of simulated particles for some cases. Also, since fewer numbers of 
particles are required for a given flow, the range of Knudsen numbers readily simulated can be extended further into 
the continuum regime. One problem encountered in the calculation of the number of tentative collisions is how to 
define the collision cell number density. For higher density simulations, the number of collisions near the axis may 
be too high, resulting in a higher mean collision rate and lower mean free path near the axis. A modified method to 
calculate the number of tentative collisions, which takes into account the varying molecular radial weights within a 
collision cell, has been proposed and shown to be independent of how the center of a collision cell is defined. 

Another important issue when performing axially symmetric simulations is the inclusion of a duplication buffer. 
A modification to the cloned molecules has been presented where the radial and axial velocities are unchanged while 
the out-of-plane velocity is reversed, thus preserving the kinetic energy of the molecule duplicated while changing 
nothing about the axially symmetric nature of the problem. 

It has also been shown that the method by which one treats the grid and molecular movements plays an important 
role in the calculation of pressure along the axis. Comparisons have been made of two traditional methods to move 
the molecules through their respective grid scheme, each having benefits and detractors. A new, hybrid grid 
/molecule translation scheme has been proposed combining the benefits of both traditional schemes. This new 
scheme uses reflection planes to enforce axial symmetry, as DAC does, but the reflection plane angle is a function of 
the collision cell radial position and height. Due to the variable reflection plane angle, the grid merit parameter has 
been reduced compared to the constant angle that DAC uses. Although the merit parameter is higher than that of the 
DS2 implementation, results of the new scheme have shown a much better distribution of pressure along the axis. 
While the modifications presented herein have alleviated the issues discussed above, further development of these 
algorithms must continue to account for the variation of molecular weights within cells. 
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